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ABSTRACT 

In  this  paper  we  review  the  development  of  the  shock-capturing  methodology, 
paying  special  attention  to  the  increasing  nonlinearity  in  its  design  and  its  relation 
to  interpolation.  It  is  well-known  that  high-order  approximations  to  a  discontinuous 
function  generate  spurious  oscillations  near  the  discontinuity  (Gibbs  phenomenon). 
Unlike  standard  finite-difference  methods  which  use  a  fixed  stencil,  modern  shock¬ 
capturing  schemes  use  an  adaptive  stencil  which  is  selected  according  to  the  local 
smoothness  of  the  solution.  Near  discontinuities  this  technique  automatically  switches 
to  one-sided  approximations,  thus  avoiding  the  use  of  discontinuous  data  which  brings 
about  spurious  oscillations. 
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1  Introduction 


In  this  paper,  we  describe  and  analyze  numerical  techniques  that  are  designed  to 
approximate  weak  solutions  of  hyperbolic  systems  of  conservation  laws  in  several 
space  dimensions.  For  sake  of  exposition,  we  shall  describe  these  methods  as  they 
apply  to  the  pure  initial  value  problem  (IVP)  for  a  one-dimensional  scalar  conservation 
law 

ut  +  f(u)  i  =  0,  u(x,0)  =  u0(x).  (1.1) 

To  further  simplify  our  presentation,  we  assume  that  the  flux  f(u )  is  a  convex  function, 
i.e.,  /"(it)  >  0  and  that  the  initial  data  uq(x)  are  piecewise  smooth  functions  which 
are  either  periodic  or  of  compact  support.  Under  these  assumptions,  no  matter  how 
smooth  Uo  is,  the  solution  u(x,t)  of  the  IVP  (1.1)  becomes  discontinuous  at  some 
finite  time  t  =  tc.  In  order  to  extend  the  solution  for  t  >  tc,  we  introduce  the  notion 
of  weak  solutions,  which  satisfy 

d  rb 

—  J  u  dx  +  f(u(b,t))  -  f(u(a,t))  =  0  (1.2a) 

for  all  6  >  a  and  t  >  0.  Relation  (1.2a)  implies  that  u(x,  t )  satisfies  the  PDE  in  (1.1) 
wherever  it  is  smooth,  and  the  Rankine-Hugoniot  jump  relation 

f{u(y  +  0,  t))  -  f{u(y  -  0,  <))  =  [u(y  +  0,  t)  -  u(y  -  0,  t)]^  (1.26) 

across  curves  x  =  y(t)  of  discontinuity. 

It  is  well-known  that  weak  solutions  are  not  uniquely  determined  by  their  initial 
data.  To  overcome  this  difficulty,  we  consider  the  IVP  (1.1)  to  be  the  vanishing 
viscosity  limit  e  J.  0  of  the  parabolic  problem 

(ue)t  +  f(ue)x  =  euexx  ue(x,0)  =u0{x),  (1.3a) 

and  identify  the  unique  “physically  relevant”  weak  solution  of  (1.1)  by 


u  =  limu*. 

ejO 


(1.36) 


The  limit  solution  (1.3)  can  be  characterized  by  an  inequality  that  the  values  ui  — 
u(y  —  0,t),uR  —  u(y  +  0,f)  and  s  =  dy/dt  have  to  satisfy;  this  inequality  is  called  an 
entropy  condition;  admissible  discontinuities  are  called  shocks.  When  f(u)  is  convex, 
this  inequality  is  equivalent  to  Lax’s  shock  condition 


o(ul)  >  3  >  a(uR) 


(1.4) 


where  a(u)  —  f'(u)  is  the  characteristic  speed  (see  [8]  for  more  details). 

We  turn  now  to  describe  finite  difference  approximations  for  the  numerical  solution 
of  the  IVP  (1.1).  Let  v"  denote  the  numerical  approximation  to  u(xj,tn)  where 
Xj  =  jh ,  tn  =  nr;  let  t/j,(x,  t)  be  a  globally  defined  numerical  approximation  associated 
with  the  discrete  values  {v”},oo  <  j  <  oo,n  >  0. 


or 


M 

□ 

□ 


n/ 
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The  classical  approach  to  the  design  of  numerical  methods  for  partial  differential 
equations  is  to  obtain  a  solvable  set  of  equations  for  {f"}  by  replacing  derivatives  in 
the  PDE  by  appropriate  discrete  approximations.  Therefore,  there  is  a  conceptual 
difficulty  in  applying  classical  methods  to  compute  solutions  which  may  become  dis¬ 
continuous.  Lax  and  Wendroff  [9]  overcame  this  difficulty  by  considering  numerical 
approximations  to  the  weak  formulation  (1.2a)  rather  than  to  the  PDE  (1.1).  For  this 
purpose,  they  have  introduced  the  notion  of  schemes  in  conservation  form: 

v-+1  =  -  A(7i+*  - 7i-j)  =  (Eh  •  lO*  (1.5a) 

here  A  =  r/h  and  fi+ 1  denotes 


/»+ 1  —  f(V?-k+U  •  •  •  iV?+k)’>  (1.56) 

/( Wi, . . .  ,w2k )  is  a  numerical  flux  function  which  is  consistent  with  the  flux  f(u),  in 
the  sense  that 

J(u,u,...u)  =  /(u);  (1.5c) 

Eh  denotes  the  numerical  solution  operator.  Lax  and  Wendroff  proved  that  if  the 
numerical  approximation  converges  boundedly  almost  everywhere  to  some  function 
u,  then  u  is  a  weak  solution  of  (1.1),  i.e.,  it  satisfies  the  weak  formulation  (1.2a). 
Consequently  discontinuities  in  the  limit  solution  automatically  satisfy  the  Rankine- 
Hugoniot  relation  (1.2b).  We  refer  to  this  methodology  as  shock-capturing  (a  phrase 
coined  by  H.  Lomax). 

In  the  following,  we  list  the  numerical  flux  function  of  various  3-point  schemes 
(k  =  1  in  (1.5b)): 

(i)  The  Lax-Friedrichs  scheme  [7] 

f(wuw2)  =  ^[f(wr)  +  f(w2)  -  j(w2  -  u>i)J  (1.6) 

(ii)  Godunov’s  scheme  [1] 

J(wuw2)  =  f(V{0]Wuw2))-,  (1.7  a) 


here  V(x/t\w1}w2 )  denotes  the  self-similar  solution  of  the  IVP  (1.1)  with  the  initial 
data 

:  tl  ■  w 


(iii)  The  Cole-Murman  scheme  [12]: 


f(w!,w2)  =  ~[f(w i)  +  f(w2)  -  \a(wl,w2)\(w2  -  toi)] 


(1.8a) 


where 


a(w  uw2)  =  <  *»r*i 


if  wx  ^  w2 


a(wi )  if  W\  =  w2 


(1.86) 
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(iv)  The  Lax-Wendroff  scheme  [9]: 

l(wi,w2)  =  ^{f(w i)  +  f(w2)  -  X a  [fM  -  f(w  1)]}.  (1.9) 

Let  E(t )  denote  the  evolution  operator  of  the  exact  solution  of  (1.1)  and  let  Eh 
denote  the  numerical  solution  operator  defined  by  the  RHS  of  (1.5a).  We  say  that  the 
numerical  scheme  is  r-th  order  accurate  (in  a  pointwise  sense)  if  its  local  truncation 
error  satisfies 

E(t)  ■  u  -  Eh  ■  u  --  0(hT+1)  (1.10) 

for  all  sufficiently  smooth  u\  here  r  =  0 {K).  If  r  >  0,  we  say  that  the  scheme  is 
consistent. 

The  schemes  of  Lax-Friedrichs  (1.6),  Godunov  (1.7),  and  Cole-Murman  (1.8)  are 
first  order  accurate;  the  scheme  of  Lax-Wendroff  (1.9)  is  second  order  accurate. 

We  remark  that  the  Lax-Wendroff  theorem  states  that  if  the  scheme  is  convergent, 
then  the  limit  solution  satisfies  the  weak  formulation  (1.2b);  however,  it  need  not  be 
the  entropy  solution  of  the  problem  (see  [4]).  It  is  easy  to  see  that  the  schemes  of 
Cole-Murman  (1.8)  and  Lax-Wendroff  (1.9)  admit  a  stationary  “expansion  shock” 
(i.e.,  }{ui)  =  /( ur )  with  a{ui)  <  a(un))  as  a  steady  solution.  This  problem  can  be 
easily  rectified  by  adding  sufficient  numerical  dissipation  to  the  scheme  (see  [11]  and 

[3])- 

2  Interpolator y  Schemes  and  Linear  Discontinu¬ 

ities 

Let  us  consider  the  constant  coefficient  case  f(u)  =  au,  a  =  const,  in  (1.1),  i.e., 

ut  +  aux  =  0,  u(x,  0)  =  it0(x),  (2.1a) 

the  solution  to  which  is 

u(x,  t)  =  u0(x  —  at).  (2.1  b) 

In  this  case  the  schemes  (1.6)  -  (1.9)  take  the  form 

v"+,=  £  =  (B*  ■  »“)i  (2  2) 

l=-K~ 

where  v  =  Aa  is  the  CFL  number.  The  coefficients  C/( v)  are  independent  of  the 
numerical  solution  vn\  this  makes  Eh  a  linear  operator. 

We  say  that  the  numerical  scheme  Eh  is  (linearly)  stable  if 

||(£h)l  <  C  for  0  <  nr  <  T,  t  =  0(h).  (2.3 a) 

In  the  constant  coefficient  case  the  scheme  is  stable  if  and  only  if  it  satisfies  von 
Neumann’s  condition 

K+ 

|  £  Ct(v)l^\  <  1  for  all  0<£<tt.  (2.36) 

t-  -K- 
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It  is  easy  to  see  that  all  the  3  point  schemes  (1.6)  -  (1.9)  are  stable  under  the  OFL 
condi  tion 

M  =  |Aa|  <  1.  (2.3c) 

The  notion  of  stability  (2.3a)  is  related  to  convergence  through  Lax’s  equivalence 
theorem,  which  states  that  a  consistent  linear  scheme  is  convergent  if  and  only  if  it 
is  stable  (see  [13]  for  more  details). 

Let  us  denote  by  S*  the  stencil  of  (r  -f  1)  successive  points  starting  with 


—  {xt>  ^t+i )  •  •  •  j  *E«+r  (2.4a) 

let  P(x;  Sf,  u)  denote  the  unique  polynomial  of  degree  r  interpolating  the  (r-fl)  values 
of  u  on  this  stencil  and  let  Q(x;  u )  denote  the  piecewise  polynomial  interpolation  of 
u 

Q(X]U)  =  P(X\ Sl^yU)  Xj_i  <  X  <  Xj  (2.46) 

We  refer  to  the  numerical  scheme 

u"+1  =  Q(xj  -  ar;vn )  (2.4c) 


as  interpolatory  scheme.  Clearly,  the  interpolator  scheme  (2.4)  is  r-th  order  a.  -.urate. 
When  Q(x;  vn)  is  the  piecewise  linear  interpolation  of  vn  (i.e.,  r  =  1  ,i(j)  —  j  —  1 
in  (2.4b))  then  (2.4c)  is  the  first-order  accurate  upwind  scheme;  in  the  constant 
coefficient  case  this  scheme  is  identical  to  those  of  Godunov  (1.7)  and  Cole-Murman 
(1.8). 


Next  let  us  assume  a  >  0  and  consider  the  second  order  case  r  —  2  in  which 
Q(x ;  u")  is  a  piecewise-parabolic  interpolation  of  vn.  There  are  two  different  choices  of 
stencil  in  (2.4):  Taking  Q  in  [xj_a,  Xj]  to  be  the  parabola  through  Sj_1  =  {xj-i,Xj,  Xj+1 } 
(i.e.,  i(j )  =  j  —  1)  results  in  the  Lax-Wendroff  scheme  (1.9);  taking  Q  in  [xj_i,Xj] 
to  be  the  parabola  through  Sj_2  =  {xj_2,  Xj_i,  Xj}  (i.e.,  ( i{j )  ~  j  -  2)  results  in  the 
second-order  upwind  scheme. 

We  turn  now  to  consider  the  application  of  these  schemes  to  the  step  function 
H{x) 


H(x) 


0x<0  _  f  0  j < 0 

1  x  >  0  ’ ~  \  1  j  >  1  ’ 


(2.5a) 


For  the  first  order  upwind  scheme  we  get  that 


Q(x-,H) 


'0  x  <  0 
<  x/h  0  <  x  <  h  ; 
1  h  <  x 


for  the  Lax-Wendroff  scheme 


Q(x-H) 


0  x  <  —h 

2l(l  +  !)  —h<x<0  . 

l-5(l-!K2-f)  0<x<h  • 

(  1  h  <  x 


(2.56) 


(2.5c) 
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for  the  second  order  upwind  scheme  we  get  that 


Q(x-H) 


0 

11(1 +  S) 

i  +  i(!-i)(2 


x  <  0 
0  <  x  <  h 
l)  h<x<2h  ‘ 
2  h  <  x 


(2.5  d) 


We  observe  that  Q  in  (2.5b)  is  a  monotone  function  of  x;  consequently  the  numerical 
solution  by  Godunov’s  scheme  to  these  data  is  also  monotone.  On  the  other  hand 
Q  for  the  second  order  schemes  (2.5c)  -  (2.5d)  is  not  a  monotone  function.  For  the 
Lax-Wendroff  scheme  Q  is  negative  in  —  h  <  x  <  0  and  has  a  minimum  of  -0.125; 
similarly  for  the  second  order  upwind  scheme  Q  is  larger  than  1  in  h  <  x  <  2h 
with  a  maximum  of  1.125.  This  observation  explains  the  Gibbs-like  phenomenon  of 
generating  spurious  oscillations  in  calculating  discontinuous  data  with  these  second 
order  schemes. 

We  say  that  the  scheme  Eh  is  monotonicity  preserving  if 


v  monotone  =>■  Eh  ■  v  monotone. 


(2.6) 


Clearly  the  numerical  solution  of  a  monotonicity  preserving  scheme  to  initial  data  of  a 
step-function  is  always  monotone  and  therefore  the  discontinuity  propagates  without 
generating  spurious  oscillations. 

Godunov  has  shown  that  the  linear  scheme  (2.2)  is  monotonicity  preserving  if  and 
only  if 

Ce(v)>  0,  -K.<i<K+]  (2.7) 

this  implies  that  a  monotonicity-preserving  scheme  which  is  linear  is  necessarily  only 
first-order  accurate.  It  took  some  time  to  realize  the  Godunov’s  monotonicity  the¬ 
orem  does  not  mean  that  there  are  no  high-order  accurate  monotonicity  preserving 
schemes;  it  only  means  that  there  are  no  such  linear  ones.  Hence  high-order  accurate 
mono  tonicity- preserving  schemes  are  nonlinear  in  an  essential  way. 

The  second-order  accurate  schemes  mentioned  above  are  linear  because  the  choice 
of  the  stencil  (2.4)  is  fixed.  Let  us  consider  now  a  piecewise-quadratic  interpolation 
which  is  made  nonlinear  by  an  adaptive  selection  of  the  stencil  in  (2.4b).  For  the 
interval  [xy_i,Xj]  let  us  consider  the  two  stencils  Sj_2  —  ixj-2>  xj-i,  xj}  and  = 
Xj+i},  and  select  the  one  in  which  the  interpolant  is  smoother.  If  we  measure 
the  smoothness  of  u  by  the  second  derivative  of  the  corresponding  parabola  we  select 


j~  2 
j  ~  1 


if  |^P(x;S?_2,u)|  <  |£P(s;5j_1>Ti)| 

otherwise 


(2.8a) 


When  we  apply  this  selection  of  stencil  to  the  step-function  H(x )  (2.5a)  we  get  that 
for  [x_i,x0]  we  choose  the  stencil  S2  2  =  {x_2,  x_lt  x0}  for  which  P(x;  Sl2,  H)  =  0;  for 
the  interval  [xi,  x2]  we  choose  the  stencil  52  =  {xi,  x2,  x3}  for  which  P(x;  S2,  H )  =  1. 
As  is  evident  from  comparing  (2.5c)  and  (2.5d)  it  does  not  matter  which  stencil  we 
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assign  to  [x0)*i]  since  both  parabolae  are  monotone  there;  with  (2.8a)  we  select 
for  [x0,  sj].  Thus  we  get  in  (2.4) 


(  0  x  <  0 


Q(x;H)={ 


if  Cl  +  s)  0  <  x  <  h 

1  h  <  x 


(2.8  b) 


which  is  a  monotone  function  of  x  although  it  is  actually  a  piecewise-quadratic  poly¬ 
nomial. 

The  use  of  an  adaptive  stencil  is  the  main  idea  behind  the  Essentially  Non- 
Oscillatory  (ENO)  schemes  to  be  described  later  in  this  paper.  It  extends  to  high 
order  of  accuracy  in  a  straightforward  manner:  For  r-th  order  accuracy  we  consider 
for  [xj_i,  Xj]  the  r  stencils  Sj_T,  Sj_r+1, . . . ,  Sj_v  We  choose  i(j)  in  (2.4b)  to  be  the 
one  which  minimizes 

\~P(x;Sri}u)\  for  i  —  j  —  r,. . .  ,j  —  1.  (2.9) 


3  Total  Variation  Stability  and  TVD  Schemes 

An  immense  body  of  work  has  been  done  to  find  out  whether  stability  of  constant 
coefficient  scheme  with  respect  to  all  “frozen  coefficients”  associated  with  the  problem, 
implies  convergence  in  the  variable  coefficient  case  and  in  the  nonlinear  case. 

In  the  variable  coefficient  case,  where  the  numerical  solution  operator  is  linear 
and  Lax’s  equivalence  theorem  holds,  it  comes  out  that  the  stability  of  the  variable 
coefficient  scheme  depends  strongly  on  the  dissipativity  of  the  constant  coefficient 
one,  i.e.,  on  the  particular  way  it  damps  the  high-frequency  components  in  the  Fourier 
representation  of  the  numerical  solution. 

In  the  nonlinear  case,  under  assumptions  of  sufficient  smoothness  of  the  PDE,  its 
solution  and  the  functional  definition  of  the  numerical  scheme,  Strang  proved  that 
linear  stability  of  the  first  variation  of  the  scheme  implies  its  convergence;  we  refer 
the  reader  to  [13]  for  more  details. 

In  the  case  of  discontinuous  solutions  of  nonlinear  problems,  linearly  stable  schemes 
are  not  necessarily  convergent;  when  such  a  scheme  fails  to  converge,  we  refer  to  this 
case  as  “nonlinear  instability.”  The  occurrence  of  a  nonlinear  instability  is  usually 
associated  with  insufficient  numerical  dissipation  which  triggers  exponential  growth 
of  the  high-frequency  components  of  the  numerical  solution. 

The  following  theorem  states  that  a  stronger  sense  of  stability,  namely  uniform 
boundedness  of  the  total  variation  of  the  numerical  solution,  does  imply  convergence 
to  a  weak  solution. 

Theorem  3.1.  Let  Vh  be  a  numerical  solution  of  a  conservative  scheme  (1.5). 

(i)  If 

TV(vh(-,t))<C-TV(u0 )  (3.1) 

where  TV(  )  denotes  the  total  variation  in  x  and  C  is  a  constant  independent  of  h 
for  0  <  t  <  T,  then  any  refinement  sequence  h  — ►  0  with  r  =  0(h)  has  a  convergent 
subsequence  hj  — >  0  that  converges  in  L]°c  to  a  weak  solution  of  (1.1). 
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(ii)  If  Vh  is  consistent  with  an  entropy  inequality  which  implies  uniqueness  of  the 
IVP  (1.1),  then  the  scheme  is  convergent  (i  e. ,  all  subsequences  have  the  same  limit, 
which  is  the  unique  entropy  solution  of  the  IVP  (1.1)). 

We  say  that  the  scheme  Eh  is  Total  Variation  Diminishing  (TVD)  if 


TV{Eh  ■  v)  <  TV(v)  (3.2) 

where 

TV{w)  =  Yl  K+i  -  wi\-  (3-3) 

j 

Clearly  TVD  schemes  satisfy  (3.1)  with  (7  =  1  and  therefore  are  TV  stable. 

In  [2]  we  have  shown  that  if  the  scheme  can  be  written  in  the  form 


v”+1  =  v?  +  <7+  x  Aj+kvn  -  C;  t  A  _ii 

J  J  J+j  "a  }~2  J  1 


(3.4a) 


where  satisfy  for  all  j 


<?**><>,  ciJ  +  e-.  <1 


j+l 


J+ 


(3.46) 


then  the  scheme  is  TVD;  here  Ai+ivn  =  -  v".  Applying  this  lemma  to  the 

general  scheme 


v 


,n+i  _ 


1 


/ j+k  ~  fj+1  ) 

we  get  that  if  A q  satisfies 

A|ai+i|  <  Agi+i  <  1 
then  the  scheme  (3.5)  is  TVD;  here 

—  /;+i  ~  fi 

Bi+*  ~  ' 


(3.5a) 

(3.56) 

(3.6a) 

(3.66) 


This  shows  that  the  Cole-Murman  scheme  (1.8)  for  which  q  =  |a|  is  TVD  subject  to 
the  CFL  restriction  A|aJ+i|  <  1. 

Using  conditions  (3.4b)  it  is  possible  to  construct  TVD  schemes  which  are  second- 
order  accurate  in  the  Lj-sense  (see  [2]  and  [14]).  However,  TVD  schemes  are  at 
most  second-order  accurate  (see  [5]).  In  order  to  design  higher-order  accurate  shock 
capturing  schemes  we  introduce  the  notion  of  Essentially  Non-Oscillatory  (ENO) 
schemes. 


4  ENO  Schemes 

In  this  section  we  describe  high-order  accurate  Godunov-type  scher'..:  which  are  a 
generalization  of  Godunov’s  scheme  (1.7)  and  van  Leer’s  MUSCL  scheme  [10]. 
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We  start  with  some  notations:  Let  {/,}  be  a  partition  of  the  real  line;  let  A(I) 
denote  the  interval-averaging  (or  “cell-averaging”)  operator 

A(J)  -w  =  ~  jw(y)dy\  (4.1) 

let  TDj  —  A(Ij)  •  w  and  denote  w  =  {ur,}.  We  denote  the  approximate  reconstruction 
of  w(x)  from  its  given  cell-averages  {w:}  by  R(x;w).  To  be  precise,  R(x-,w)  is  a 


piecewise-polynomial  function  of  degree  (r  —  1),  which  satisfies 

(i)  R(x\w)  —  w(x)  +  0(hT)  wherever  w  is  smooth  (4.2a) 

(ii)  A(/j)  '  E.(-\w)  =  Wj  (conservation).  (4-26) 

Finally,  we  define  Godunov-type  schemes  by 

u"+1  =  A(lj)  ■  E(t)  •  R(-.vn)  =  (Eh  ■  vn),  (4.3a) 

v°  =  A(Ij)u0 ;  (4.36) 

here  E{t)  is  the  evolution  operator  of  (1.1). 


In  the  scalar  case,  both  the  cell-averaging  operator  A{I})  and  the  solution  operator 
E(t)  are  order-preserving,  and  consequently  also  total-variation  diminishing  (TVD); 
hence 

TV (Eh  ■  w)  <  TV(R(-tW)).  (4.4) 

This  shows  that  the  total  variation  of  the  numerical  solution  of  Godunov-type 
schemes  is  dominated  by  that  of  the  reconstruction  step. 

We  turn  now  to  describe  the  recently  developed  essentially  non-oscillatory  (ENO) 
schemes  of  [5,  6],  which  can  be  made  accurate  to  any  finite  order  r.  These  are 
Godunov-type  schemes  (4.3)  in  which  the  reconstruction  R(x\w),  in  addition  to  re¬ 
lations  (4.2),  also  satisfies 

TV{R(-w))  <  TV(w)  +  0{h 1+p),  p  >  0  (4.5) 

for  any  piecewise-smooth  function  w(x),  Such  a  reconstruction  is  essentially  non- 
oscillatory  in  the  sense  that  it  may  not  have  a  Gibbs-like  phenomenon  at  jump- 
discontinuities  of  w(x),  which  involves  the  generation  of  0(1)  spurious  oscillations 
(that  are  proportional  to  the  size  of  the  jump);  it  can,  however,  have  small  spurious 
oscillations  which  are  produced  in  the  smooth  part  of  w(x),  and  are  usually  of  the 
size  0(hT)  of  the  reconstruction  error  (4.2a). 

When  we  use  an  essentially  non-oscillatory  reconstruction  in  a  Godunov-type 
scheme,  it  follows  form  (4.4)  and  (4.5)  that  the  resulting  scheme  (4.3)  is  likewise 
essentially  non-oscillatory  (ENO)  in  the  sense  that  for  all  piecewise-smooth  function 
w(x) 

TV(Eh  ■  w)  <  TV(w)  +  0(h1+p),  p>  0;  (4.6) 

i.e.,  it  is  “almost  TVD.”  Property  (4.6)  makes  it  reasonable  to  believe  that  the  total 
variation  of  the  numerical  solution  is  uniformly  bounded.  We  recall  that  by  Theorem 
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3.1,  this  would  imply  that  the  scheme  is  convergent  (at  least  in  the  sense  of  having 
convergent  subsequences).  This  hope  is  supported  by  a  very  large  number  of  numerical 
experiments. 

Next  we  describe  one  of  the  techniques  to  obtain  an  ENO  reconstruction.  To 
simplify  our  presentation  we  assume  that  {Ij}  is  a  uniform  partition 


Ij  —  (xj_  Xj  —  jh. 


Given  cell  averages  {w3}  of  piecewise-smooth  function  w(x ),  we  observe  that 


hwj  =  f  w(y)dy  —  W(xj)  —  W(xj_i)  (4.7a) 

Jxj-l 

where  x 

W(x)  =  f  u >(y)dy  (4-76) 

Jx  0 

is  the  primitive  function  of  w(x).  Hence  we  can  easily  compute  the  point  values 
{W(xj))  by  summation 

t 

W(xt)  =  h  ^2  Wj-  (4.7c) 

J  =  «o 

Once  we  have  computed  the  point  values  of  the  primitive  function  we  use  the  ENO  in¬ 
terpolation  technique  (2.4),  (2.9)  to  obtain  Q(x\W),  an  r-th  order  piecewise- 
polynomial  interpolation  of  W ,  i.e., 

Q(x;  IV)  =  P(x;  S^,  W)  for  x}_\  <  x  <  Xj  (4.8a) 


where  P(x;5lr,  W)  is  the  unique  r-th  degree  polynomial  which  interpolates  W  over 
the  stencil  5,r  =  {x^,  x,+i, . . . ,  x,+r},  and  i(j)  is  chosen  so  that 


(4.86) 


We  define  R(x,w)  by 

R(x-W)  =  ^Q(x;W).  (4.9) 

We  observe  that  if  w{x)  is  smooth  in  (x;_ i,x;)  then  for  h,  sufficiently  small  the 
algorithm  (4.8b)  will  select  a  stencil  in  which  w(x)  is  smooth.  It  follows  then 
from  standard  interpolation  theorems  that 

R(z;w)  =  ^P(x;  ^(j)>  W)\  =  +  0(hr)  =  w(x)  +  0(hr)  (4.10) 

which  is  property  (4.2a).  Furthermore  (4.10)  holds  in  every  interval  except  for  those 
in  which  w(x)  has  a  discontinuity.  As  we  have  seen  in  the  examples  (2.5)  and  (2.8b) 
the  Gibbs-phenomenon  is  associated  with  intervals  near  the  discontinuity  and  not 
with  the  interval  that  contains  the  discontinuity.  This  is  why  the  reconstruction  (4.8) 
-  (4.9)  satisfies  the  ENO  property  (4.5);  in  [2]  we  show  that  the  second-order  accurate 
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ENO  scheme  is  actually  TVD.  The  conservation  property  (4.2b')  follows  directly  from 
the  definition  (4.9): 


A^Ri-w) 


\J’  w^dx  =  w)  -  Qix3- 1;  w)] 

^-[VE(xj)  -  W{x]^l)}  =  wr 


(4.11) 


The  abstract  scheme  (4.3)  can  be  written  in  the  standard  conservation  form  (1-5). 
To  do  so  let  us  denote  by  v(x,t)  the  solution  in  the  small  of  the  IVP 


(u)t  +  f(v)x  =  0 
o(x,  0)  —  R{x\  vn ) 


0  <  t  <  T 


(4-12) 


and  integrate  this  PDE  over  /_,  x  [0,  r] ;  using  the  divergence  theorem  and  (4.2b)  wc 
get  that  un+1  in  (4.3)  can  be  expressed  by 


,«+ 1 


(4.13a) 


where 

fj+i  =  ~  [  f{v(xj,t))dt  (4.136) 

3  T  JO 

In  the  first-order  case  the  scheme  (4.13)  is  identical  to  Godunov’s  scheme  and  the 
numerical  flux  (4.13b)  can  be  expressed  in  a  closed  form  by  (1.7b).  For  higher  order 
schemes  we  use  a  numerical  flux  which  is  an  appropriate  approximation  to  (4.13b) 
(see  [6]  for  more  details). 

We  remark  that  the  ENO  schemes  are  related  to  the  interpolatory  schemes  of  Sect. 
2  as  follows:  In  the  constant  coefficient  case  a  fixed  choice  of  stencil  (i.e.,  i(j)  —  j 
=  constant  in  (4.8a))  results  in  the  interpolatory  scheme  (2.4)  corresponding  to  the 
same  choice  of  stencil. 
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